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PHASE RETRIEVAL IN PHASE CONTRAST IMAGING 

10 Field of the Invention 

This invention relates generally to the observation of structural features of objects 
utilising penetrating radiation such as x-rays. In particular, the invention relates to the 
derivation of images of the phase change introduced by an object in penetrating radiation 
incident on the object, from a two-dimensional intensity record of the penetrating 
15 radiation after it has traversed the object. The invention may be extended to retrieve 
separate phase and absorption data from a set of radiographic measurements. 

Background Art 

The present applicant's international patent publications WO 95/05725 
20 (PCT/AU94/00480) and WO 96/31098 (PCT/AU96/00178) disclose various 
configurations and conditions suitable for differential phase-contrast imaging using hard 
x-rays. Other earlier disclosures of interest are to be found in Soviet patent 1402871 and 
in US patent 5319694. Differential phase-contrast imaging shows great promise for 
viewing the internal structure of objects for which traditional absorption-contrast 
25 radiography is of limited or no value because of very weak absorption contrast. This is 
the case, for example, with soft tissue within the human body. 

The practical issue of optimally and efficiently deriving the phase-contrast image 
for an object from the actual record at the detector is addressed in two related papers by 
30 Nugent et al Phys. Rev. Lett. 77, 2961-2964 (1996); J. Opt. Soc. Am.. A13, 1670-82 
(1996) and references therein. In these papers, it has been demonstrated that with 
monochromatic plane-wave x-radiation as in the configurations of WO 95/05725 and US 
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5319694, the retrieval of phase information from measurements of the propagation of 
intensity can be based on treating the propagation of the modified radiation field whose 
characteristics reflect the phase modifying effects of the object. A two-dimensional 
recording of the intensity of the penetrating radiation after it has traversed the object is 
5 the result of variations in the local direction of propagation of the radiation arising from 
variations in local refractive index, typically an indication of a boundary or rapid 
variation in electron density within the object or of a thickness variation. The 
aforementioned articles by Nugent et al utilise a treatment of the propagation of a plane 
monochromatic electromagnetic wave based on Maxwell's equations to derive a 

10 transport-of-intensity equation and propose solutions of this equation to derive a phase- 
contrast image from the intensity record. These suggested solutions to the transport-of- 
intensity equation involve expanding the phase in orthogonal functions. The kind of 
function chosen depends on the shape of the sample, and thus Zernike polynomials are 
adequate for a circular shape whilst a Fourier expansion is most suitable for a square- 

15 shaped sample. 

The aforementioned international patent publication WO 96/31098 discloses an 
in-line phase-contrast imaging configuration utilising a substantially point source and a 
two-dimensional x-ray imaging detector spaced from the object. It is demonstrated in the 

20 application that, in contrast to previous phase-contrast imaging configurations, a point 
source may be utilised, and moreover that the source may be broadly polychromatic 
provided its radiation has high lateral spatial coherence, which in practical terms 
indicates a maximum source diameter (s) dependent upon the source to object distance 
(R,). The larger the source-object distance or the smaller the source size, the greater the 

25 lateral spatial coherence (see Wilkins et al Nature 384 335-8 (1996). A consequence of 
these disclosures in WO 96/31098 is that the approach proposed is more closely related 
to traditional methods used for absorption-contrast radiography and should be easier to 
implement than earlier proposals. This method of phase-contrast imaging is especially 
advantageous in the hard x-ray region where the lack of suitable lens elements make 

30 other techniques conventionally used in visible light and soft x-ray microscopy 
unsuitable. 
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It is an object of the present invention, at least in one or more embodiments, to 
provide a method of obtaining a phase-contrast image from a two-dimensional intensity 
record where the penetrating radiation substantially emanates from a point-like source. 
In one or more embodiments, it is a particular objective to provide a method adaptable to 
5 the extraction of phase and absorption-contrast information from radiographic images 
recorded with a microfocus source which need not be highly monochromatic. 



Summary of the Invention 

In certain aspects, the invention involves the concept of obtaining two or more 
10 intensity records at a common finite distance after the radiation has emerged from the 
object, for respective different energy distributions of the radiation. 



In one or more other aspects, the invention entails an appreciation that, with 
certain features, a point-like source configuration also lends itself to an approach based 
15 on the solution of a differential transport-of-intensity equation, albeit a different one 
from that used by others for the plane-wave case. 

The invention is especially useful for separating out and retrieving phase 
information from a typical intensity record (obtained with a microfocus radiation source) 
20 which has both phase-contrast and absorption-contrast content. 



In a first aspect, the invention provides method of obtaining an image of the 
phase change introduced by an object in penetrating radiation incident on the object, 
including: 

25 irradiating the object with penetrating radiation having high lateral spatial 

coherence; 

receiving at least a portion of said radiation at detector means after the 
radiation has emerged from the object and thereby obtaining and storing at least 
two intensity records for the received radiation each including intensity values at 
30 predetermined intervals; and 

utilising these values to derive a grid of values defining an image of the 
phase change introduced by the object in the penetrating radiation; 
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wherein said intensity records are obtained at a uniform finite distance 
after the radiation has emerged from the object, and are for respective different 
energy distributions of the detected radiation. 



5 In its first aspect, the invention further provides apparatus for obtaining an image 

of the phase change introduced by an object in penetrating radiation incident on the 
object, including: 

means to provide a source for irradiating an object with penetrating 
radiation having high lateral spatial coherence; and 

10 detector means for receiving at least a portion of said radiation after the 

radiation has emerged from the object whereby to generate at least two intensity 
records for the received radiation each including intensity values at 
predetermined intervals; 

wherein said detector means is arranged for obtaining said intensity 

15 records at a uniform finite distance after the radiation has emerged from the 

object, and energy characterising means is provided whereby said intensity 
records are for respective different energy distributions of the detected radiation. 



In one embodiment, the respective different energy distributions are obtained by 
20 altering the energy spectrum of the radiation irradiating the object. This might be 
achieved, for example, by modifying the output of the radiation source, or by pre-filter 
means. In another embodiment, the respective different energy distributions are obtained 
by providing for the detector means to provide intensity as a function of energy in a 
certain energy band or band(s). For this purpose, a two-dimensional detector means may 
25 be variably wavelength sensitive, or be preceded by a variable filter shutter means. 
Further alternatively, the image intensity may be recorded as a function of photon energy 
for each pixel over a number of ranges of x-ray energy. Advantageously, for enhanced 
resolution, multiple intensity records may be obtained for respective multiple energy 
distributions of the radiation. 

30 

In the simplest case, each energy distribution may be a particular wavelength or 
photon energy level. 
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The aforesaid derivation may include solving one or more differential transport- 
of-intensity equations relating the phase at a plane of the object to the evolution of the 
intensity distribution along the direction of propagation, utilising predetermined uniform 
5 boundary conditions. An alternative derivation includes solving Fourier-optics 
equations. Others are of course possible for particular configurations or circumstances. 

Preferably, the intensity values also reflect absorption contrast in the object, and 
the method further includes utilising these values to derive a grid of values defining an 
10 effective pure absorption-contrast image of the object. 

With the apparatus of the first aspect of the invention, there is preferably further 
included a computer program product having a set of machine readable instructions 
which, when installed in a computer having a suitable operating system and memory 
15 means, configures the computer to be operable to utilise said values to derive a grid of 
values defining an image of the phase change introduced by the object in the penetrating 
radiation. 

In a second aspect, the invention provides a method of obtaining an image of the 
20 phase change introduced by an object in penetrating radiation incident on the object, 
from one or more two-dimensional intensity records of penetrating radiation after it has 
traversed the object, the radiation being of high lateral spatial coherence when incident 
on the object and the or each record being obtained at a finite distance after the radiation 
has emerged from the object incorporating phase-perturbed components within a 
25 surrounding field of radiation either uniformly phase-perturbed or not phase-perturbed, 
the method including: 

storing intensity values from the or each record, at predetermined 
intervals; 

utilising these values and any predetermined uniform boundary conditions 
30 to derive a grid of values defining an image of the phase change introduced by 

the object in the penetrating radiation, by solving a differential transport-of- 
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intensity equation relating the phase at an exit plane at the object to the evolution 
of the intensity distribution along the direction of propagation. 

The invention also provides, in its second aspect, a method of obtaining an image 
5 of the phase change introduced by an object in penetrating radiation incident on the 
object, including: 

irradiating the object with penetrating radiation having high lateral spatial 
coherence; 

receiving at least a portion of the said radiation at a detector at one or 
10 more finite distances after the radiation has emerged from the object 

incorporating phase-perturbed components within a surrounding field of radiation 
either uniformly phase-perturbed or not phase-perturbed, and thereby obtaining 
and storing intensity values for the received radiation at predetermined intervals; 
and 

15 utilising these values and any predetermined uniform boundary conditions 

to derive a grid of values defining an image of the phase change introduced by 
the object in the penetrating radiation, by solving a differential transport-of- 
intensity equation relating the phase at an exit plane at the object to the evolution 
of the intensity distribution along the direction of propagation. 

20 

Preferably, the same intensity values are further utilised to also derive values 
defining an effective pure absorption-contrast image for the object. 

The invention further provides, in its second aspect, an apparatus for obtaining an 
25 image of the phase change introduced by an object in penetrating radiation incident on 
the object, from one or more two-dimensional intensity records of penetrating radiation 
after it has traversed the object, the radiation being of high lateral spatial coherence when 
incident on the object and the or each record being obtained at a finite distance after the 
radiation has emerged from the object incorporating phase-perturbed components within 
30 a surrounding field of radiation either uniformly phase-perturbed or not phase-perturbed, 
the apparatus including: 
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(a) means for storing intensity values from the or each record, at pre-determined 
intervals; and 

(b) means, preferably including a stored program of machine readable instructions, 
for utilising said values and any predetermined uniform boundary conditions to 

5 derive a grid of values defining an image of the phase change introduced by the 

object in the penetrating radiation a phase-contrast image for the object, by 
solving a differential transport-of-intensity equation relating the phase at a plane 
at the object, e.g. an exit plane, to the evolution of the intensity distribution along 
the direction of propagation. 

10 

The invention also provides a set of machine readable instructions which, when 
installed in computer means also having suitable operating system software and memory 
means, configures the computer means for operation as apparatus according to the 
preceding paragraph. The invention still further provides a storage medium, e.g. a 
15 magnetic disk, a CD-ROM or optical storage disk, or an internet server, in which is 
stored said set of machine readable instructions. 



The invention still further provides, in its second aspect, apparatus for obtaining 
an image of the phase change introduced by an object in penetrating radiation incident on 
20 the object, including: 

a source for irradiating an object with penetrating radiation having high 
lateral spatial coherence; 

a detector for receiving at least a portion of said radiation a finite distance 
after the radiation has emerged from the object incorporating phase-perturbed 
25 components within a surrounding field of radiation not phase-perturbed or 

uniformly phase-perturbed , whereby to generate intensity values for the received 
radiation at pre-determined intervals; and 

computer means, including a stored program of machine readable 
instructions, operable to utilise said values and any predetermined uniform 
30 boundary conditions to derive a grid of values defining an image of the phase 

change introduced by the object in the penetrating radiation, by solving a 
differential transport-of-intensity equation relating the phase at an exit plane at 
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the object to the evolution of the intensity distribution along the direction of 
propagation. 

Preferably, the computer means is further operable to utilise the same intensity 
5 values to also derive values defining an effective pure absorption-contrast image for the 
object. 

In one embodiment of the method of the second aspect of the invention, two of 
the intensity records are obtained at different distances after the radiation has emerged 

10 from the object, e.g. in two distinct image planes. In another embodiment two or more 
of the intensity records are obtained at a uniform distance, e.g. in a single image plane, 
for different energy distributions of the incident radiation on the object. In a particular 
form of the latter embodiment, one or more of said intensity records are recorded such 
that image intensity as a function of photon energy is recorded for each pixel over a 

15 number of ranges of x-ray energy. 

The transport-of-intensity equation relevant to the point source case may be 
equation (16), or in an alternative form, equation (18), hereunder, and the solution may 
be by a perturbation method. Alternatively, the equation may be solved numerically, 
20 especially when the last two terms of equation (16) have similar magnitudes. 

As employed herein, the term "penetrating radiation" includes x-rays and 
neutrons, although the invention is especially useful for x-ray radiation, and may be 
substantially monochromatic but is more typically broadly polychromatic. An 

25 application of particular interest is the range 0.5keV to 1 MeV, e.g. the hard x-ray range 
IkeV to 1 MeV. The phase perturbation by the object may be thought of as a refraction 
effect or may more rigorously be viewed as a Fresnel diffraction effect. For an object of 
finite thickness within a surrounding medium, e.g. of different refractive index, the phase 
perturbation will also be dependent on the thickness of the object in the direction of the 

30 localised wave vector. 
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The object may, e.g., be a boundary, typically a boundary exhibiting a sharp 
variation in refractive index with respect to the penetrating radiation. The invention is 
especially useful where there is weak or negligible absorption contrast for the radiation 
between the intensities of the radiation passing through respective sides of the boundary, 
5 but may generally also be applied where there is significant absorption contrast at the 
boundary. 

The boundary conditions typically do not need to be measured and may, e.g., 
include uniform Dirichlet, Neumann or periodic boundary conditions. The boundary 
10 conditions are selected to achieve a unique solution of the equation for phase, at least up 
to an arbitrary constant component. 

Preferably, the solution further utilises one or more optical conditions. Such 
conditions may include e.g. a small wavefront curvature for the incident radiation, 
15 absence of focal points between the object and image, and uniform illumination of the 
object. 

The incident penetrating radiation is not restricted to being monochromatic. For 
polychromatic radiation, the equation includes a spectrally weighted term or factor 
20 dependent on the square of the respective wavelength components. 

Brief Description of the Figures 

The invention will now be further described, by way of example only, with 
reference to the accompanying figures, in which: 
25 Figure 1 is a diagram of an x-ray optics configuration according to an 

embodiment of the invention in which two intensity records are made at different 
detection planes; 

Figure 2 is a diagram of a related co-ordinate system for the mathematical 
discussion which follows; 
30 Figure 3 is a diagram of an alternative x-ray optics configuration in which two 

intensity records are made at a common detection plane for respective different energy 
distributions of the radiation; 
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Figures 4 to 9 illustrate successive intensity records in a mathematical test case of 
a method utilising one approach to the invention; and 

Figures 10 to 15 are photomicrographs depicting a test of an alternative approach. 

5 Preferred Embodiments 

The configuration illustrated in Figure 1 includes a microfocus source S of high 
lateral spatial coherence and an x-ray two-dimensional imaging detector D, for example 
film, photostimulable phosphor plates (eg. Fuji image plates), or a two-dimensional 
electronic detector such as a charge-coupled device (CCD) array. 

10 

The term "lateral spatial coherence" herein refers to the correlation of the 
complex amplitudes of waves between different points transverse to the direction of 
propagation of the waves. Lateral spatial coherence is said to occur when each point on 
a wavefront has a direction of propagation which does not change over time. In practice, 

15 high lateral spatial coherence may, for example, be achieved by using a source of small 
effective size or by observing the beam at a large distance from the source. Generally, 
lateral coherence length d± = A,R/s, where X is the x-ray wavelength, R t is the source to 
object distance, and s is the maximum source diameter. For example, for 20 keV x-rays 
and a source to object distance of 200 mm, a source size of approximately 20(im 

20 diameter or less would typically be appropriate. The smaller the source size the better 
for the purposes of this invention, provided total flux from the source is sufficient. 
Lateral spatial coherence may need to be preserved by careful selection of the x-ray 
window of the source, e.g. such that it is of highly uniform thickness and homogeneity. 
The effects of partial spatial and temporal coherence on image contrast and resolution are 

25 described in Pogany et al [Rev. Sci. Instrum. 68 2774-82 (1997)]. 

Regions of refractive index variation transverse to the local direction of 
propagation, or thickness or density variations in the direction of propagation, can lead to 
a significant change in the local direction of propagation of the wave front passing 
30 through those regions. 
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Thus, referring again to Figure 1, a spherical wavefront Wl emanating from the 
point source S becomes distorted to W2 on passing through the object O. By recording 
the intensity of the wavefront at a sufficient distance from the sample, intensity 
variations due to rapid refractive index and thickness or density variations in the sample 
5 may be detected and their location recorded in an image. This corresponds to a form of 
differential phase-contrast imaging. The location of the imaging detector is chosen such 
that the spatial resolution of the detector is sufficient to resolve the intensity differences 
arising from the severe distortions of the wavefront and to optimise contrast subject to 
practical considerations. The values recorded at the detector also reflect absorption in 
10 the object and so the intensity values recorded at the detector also contain extractable 
absorption-contrast information. 

For reasons to be explained subsequently, it is noted here that the object is 
disposed so that radiation which emerges from the object incorporates phase-perturbed 
15 components within a surrounding field of radiation not phase-perturbed or uniformly 
phase-perturbed. 

Typically, the sharp gradients in refractive index or thickness will be imaged as 
sharp losses or rapid variations in intensity at corresponding points in the image. This 
20 feature of intensity loss or rapid variation at a given point in the image is to a first 
approximation essentially independent of wavelength and can therefore lead to very 
sharp contrast variations in the image even when a polychromatic source is used. 

This configuration has the feature that for a circular source distribution, the 
25 spatial resolution in the image is the same for all directions and is essentially determined 
by the source size. It also has the advantage that considerable magnification of the 
image is possible and so recording media such as photostimulable phosphor image plates 
may be used which have many desirable properties such as wide dynamic range and high 
sensitivity but not high spatial resolution. 

30 
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Detector D is typically coupled to suitable computing means M such as a personal 
computer utilising a 486 CPU, operating at 66 Mhz, and provided with appropriate 
memory and applications software. 

5 Computer means M stores in memory a set of intensity values from the intensity 

record at detector D, at predetermined intervals in two-dimensions. In this simple case, 
the detector is a flat planar detector and the intervals are on a uniform square grid. The 
manner in which these values are derived will depend upon the kind of detector: for 
example, the detector may be of a pixellated structure in which each pixel is sampled in a 
10 two-dimensional scan and the values fed serially to the computer memory. Alternatively, 
the record may be made and stored in the detector and scanned or sampled as required by 
the computer. 

Computer M further includes a control program by which the stored intensity 
1 5 values from the detector record and any predetermined uniform boundary conditions are 
utilised to derive a further grid of values defining an effective phase-contrast image in a 
selected plane of the object, by solving a differential transport-of-intensity equation. 

For simplicity of formal development, we now discuss the theoretical aspects for 
20 the case of a monochromatic point source, with a view to defining a simple differential 
equation, and solution, applicable to a weak phase object, as described by conditions (4) 
-(6). The case of a broadly polychromatic source will be considered subsequently. 

With reference to the co-ordinate system indicated in Figure 2, consider a 
25 monochromatic point source with wavelength X = 2k/ k located at point z = -R P R, » X> 
on the optical axis Z. The source illuminates an object which is located between the 
source and the plane (x, y, 0) orthogonal to the optical axis at point z = 0, the object 
being bounded on the right hand side by the plane z = 0 (Figure 2). The object is 
assumed to be thin, i.e. its z-dimension is much smaller than R r 

30 

Let 
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exp {ik^Jx 2 + y 2 +R] } 
^l + (x 2 +y 2 )/ R 



u 0 (x 9 y) = — L ^ o 7 == ^j a(x,y)exp{iky/(x,y)} (1) 



be the scalar complex amplitude at the object plane z = 0, where §(x f y) = k\\f(x,y) is the 
phase "aberration" appearing due to the transmission of the wave through the object and 
5 I(x>y) = a{x y y) is the corresponding intensity. 

The propagation of the amplitude u(x,y,0) = u 0 (x,y) into the source-free half-space 
z>0 is governed by the Helmholtz equation 

10 (V 2 + £ 2 )w = 0. (2) 

A solution to the Helmholtz equation, such that u(x,y,0) = u 0 {x,y), is given by the 
first Rayleigh-Sommerfeld integral. This integral can be written for z »X as 

15 u(x< ,y .z) s -JLjj QXp{ikr} *u 0 (jc, y)dxdy. (3) 
In r r 



where r = ^{x'-xf +(y*-y) 2 + z 2 . 

In order to simply the analysis, we assume that our object is small compared to its 
distance from the source and is weakly scattering, i.e. 

20 

\]f(x,y) s 0 and a(x y y) = 1 , when (x 2 + y 2 ) > d 2 , d«R l ; (4) 
I 3 : d " V(*oOl «(R , ) , ' (n+m) , m + n 2 > 0, (5) 
25 I a(x,y)\ >C > 0 and I d m x d n y a (x,y)\ « Gfcl 3 ™ d n y \\f(x,y)\ , m + n 2 > 0 (6) 
where R' = R,R 2 /(R 1 + R 2 ) and z = R 2 is the position of the "image" plane, R 2 » A,. 
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Condition (4) requires the wavefront incident on the object to have a small 
curvature. Condition (5) also has a simple physical meaning: it ensures, in particular, 
that there are no focal points between the object plane z = 0 and the image plane z = R 2 
(as I V 2 \|/| «1/R'). In the absence of condition (5) phase distribution in the image plane 
5 may relate to intensity distribution in the object plane in a more complex manner due to 
possible focal phase shifts. Condition (6) essentially requires the incident beam to 
illuminate the object uniformly and the variations of the imaginary part of the refractive 
index inside the object to be much weaker than the variations of its real part at a given x- 
ray wavelength A,. The last requirement applies both for "phase" objects, and also for 
10 mixed phase/absorption objects. 

Conditions (4 to 6) guarantee that the diffraction angles are small. Therefore, 
under these conditions the image can be calculated using the Fresnel integral: 

15 w( jc' , y , R 2 ) = - — eXp{ ^ } J J exp {ikS(x 9 y) }a( *, y)dxdy, (7) 

2k R 2 

where R = R y + R 2 and 

r/ , (jc'-jc) 2 4-(y'-y) 2 jc 2 +v 2 , , 
(X ' y) = 2R + 2 R + ¥(Xy ^ 

20 

Applying the stationary phase formula to integral (7) we obtain: 

r- V ex P {i[kR + kS(x < >1± ) + (7r 1 2) s § n j ] ? x im ^ 
u(x ,y ,/? 2 ) = 2^ / a(x ? ,yJ + <9(/; ), 

R 2 jdetS s 

25 where the sum is over all stationary points (x ,y s ) corresponding to image point (x\ y'), 
sgn S" s is the number of negative eigenvalues of the matrix S M of second-order partial 
derivatives of the phase function S evaluated at point (x s ,y s ) and det S" s is its 



(8) 



(9) 
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determinant, I 0(k' { )\ < const / k. Conditions (6) ensure that the residual term 0(k' x ) is 
much smaller than the leading terms in (9). 

Calculating partial derivatives of the phase function (8) and equating them to 
5 zero, we obtain the following equations for the stationary points {x s , y s ): 

y=M^ + R 2V ; (x s ,y s ) (10) 
where M = (R, + R 2 ) / R, 

10 

and where we used notation of the type \j/ q ' for partial derivatives : 4 / q ' = 34^/ dq. 
Equations (10) describe ray trajectories originating on the object plane at points (x , y s ) 
and ending up at a given point (jt', y', R 2 ) of the image; constant M is the coefficient of 
image magnification. Condition (5) ensures that equations (10) have only one solution, 
15 (x , y s ), i.e. there is only one ray connecting each point of the image with the 
corresponding unique point of the object. Consequently, the summation sign in (10) can 
be omitted. 

It is straightforward to calculate the matrix of second-order partial derivatives of 
20 the phase function: 



(l//?) + v;j 



(11) 



This matrix describes the variation of the curvature of the wavefront. In view of 
25 condition (5) this matrix is non-degenerate, both its eigenvalues are positive (hence, sgn 
S" = 0) and its determinant can be approximated by 

det S"=(l + R'VVVR' 2 . Rl V>l «1. (12) 
30 Now we can simplify the stationary phase formula (9): 
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u(x ,y\R 2 ) = a(x ^ y *\ exp {iS(x, ,y,)}. (13) 



5 



10 



15 



25 



The corresponding intensity distribution on the image plane is 
I(x'y } R 2 ) = M 2 I(x s , yJH-RV 2 v(x s ,y)l (14) 

Let us now use the stationary-point equations (10) to further simply equation 

(14): 

Kx s> y s ) = KM 'V - M~ l R 2 \\r x , M"V - M 'R 2 y/' y ) 

= /(ArV, M''y') - R'(V/-V\|/)(Af V, A#V) . (15) 

Substituting (15) into (14) and denoting jc = M' 1 x\ y = M' 1 y\ we obtain 

M 2 1(Mx, My } /?J-/te^[UR'VVr^R'Vlog I(x.y) Vyfx.y)]. (16) 



This equation describes the image of a weak phase object (as described by 
conditions (4)-(6)) as a function of the source-to-object and object-to-image distances. It 
20 shows that in the case of negligible absorption the image contrast is proportional to 
"focal" distance R' and to the Laplacian of the phase shifts introduced into the incident 
beam on passing through the object Note that formula (16) in the case of V\|/ = 0 and 
V\j/ = 0 trivially gives the exact value of intensity in the absence of any objects between 
the source and the image plane. 



For x-ray wavelengths away from absorption edges of the material of the object, 
the linear attenuation coefficient jll is known to be approximately proportional to X 3 and 



A 2 r 0 

¥(^y) = ~r^ i p(x,y,z)dz, (17) 
2k -°° 
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where r is the classical electron radius and p is the electron density of the object. 
These facts and equation (16) imply that in the accepted approximation the image 
contrast has a simple quadratic dependence on the wavelength. Therefore, 
polychromaticity of the source is not necessarily an obstacle to the application of this 
method and can be explicitly taken into account : the factor X 2 is then replaced by a 
spectrally weighted sum dependent on the square of the respective wavelength 
components. 

Equation (16) can be used for explicit retrieval of the phase from intensity 
measurements in two image planes, indicated at D, and D 2 in Figure 1. To see how this 
may be achieved, it is convenient to rewrite the equation in the following form: 

-V- [I(x,y) Vvfotf =F(x,y), (18) 

where F(x,y) = [M 2 I (Mx, My, R 2 ) - I(x,y)] I R'. Equation (18) relates the transverse 
derivatives of the phase on the object plane to the intensity distributions on the object 
and image planes. 

Equation (18) is similar to the plane-wave transport-of-intensity equation (TIE) 
which can be formally obtained from (18) when R 2 — > 0. There are, however, some 
important differences. First, unlike the TIE, equation (18) does not assume the distance 
between the two planes D I? D 2 where intensity is measured, to be infinitesimally small. 
Second, equation (18) takes image magnification explicitly into account, and, therefore, 
is more appropriate for phase imaging with a point source. It can be readily seen that an 
attempt to directly apply the plane-wave TIE to imaging with a spherical incident wave 
leads to a different function on the right-hand side of equation (18), namely the distance 
R' in the denominator of F(x,y) is replaced by R 2 . The resulting error is due to the extra 
wavefront curvature produced by the incident spherical wave and is negligible only when 
R 2 « R r Finally, our method of derivation of equation (18), which was based on the 
application of the stationary phase formula to the Fresnel integral, allowed us to 
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formulate the validity conditions (4)-(6) in terms of the amplitude distribution u 0 (x,y) on 
the object plane. These are effectively conditions on the optical properties of the source, 
the sample and the geometrical parameters of the imaging layout. In contrast to this, the 
conventional method for deriving the plane-wave TIE from the Helmholtz equation 
5 involves requirements on the diffracted amplitude u(x,y,z), z > 0, which are difficult to 
verify in practice. 

In order to find a unique solution to (18) one needs to define some boundary 
conditions for it. Such boundary conditions can be readily obtained using condition (4). 
10 If Q is a simply connected domain in the plane (x,y,0) containing the circle of radius d 
with the centre at the origin of coordinates, then on the boundary r of £1 we have 

VUy)lr = 0. (19) 

15 In fact, one can define not only the Dirichlet boundary conditions (19), but any 

uniform boundary conditions on r, e.g. the Neumann or the periodic ones. This situation 
is very different from that encountered in astronomical adaptive optics, where the plane- 
wave TIE is also used for phase retrieval. There the distorted wavefront is not 
surrounded by the primary unperturbed or uniformly perturbed one as in the present 

20 situation and, therefore, boundary conditions must always be measured directly. 

In the case where the primary unperturbed or uniformly perturbed wavefront is 
not measured, a procedure of artificially blurring the image at the edges of the image 
field or of iteratively retrieving the phase and recalculating the intensity distribution in 
25 the region around the edges of the object to effectively create uniform boundary 
conditions might be used. 

Equation (18) with boundary conditions (19) has a unique solution for phase, 
provided I(x,y )>C 2 >0 in £2, where C is a constant defined in (6), while solutions to 
30 similar problems with the Neumann and periodic boundary conditions are unique up to 
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an arbitrary additive constant The requirement for the intensity to be non-zero is 
satisfied due to conditions (6). 

Solutions to the problem (18)-(19) and similar problems with the Neumann and 
5 periodic boundary conditions are stable with respect to errors in the right-hand side 
function F(x,y). Note, however, that errors in M 2 I(Mx, My, R 2 ) and I(x,y) should 
preferably be small compared not only to these intensities, but to their difference which 
has to be much smaller than the intensities themselves according to equation (16) and 
condition (5). Therefore, it is preferred, strongly, that the measured intensity data be 
10 fairly "clean" from noise to optimise the reliability of quantitative retrieval of the phase. 

In the case of uniform intensity, I(x,y) = I 0 on the object plane a simple solution 
to equation (18) with the periodic boundary conditions in the square (~D,D) x (-D,D), D 
> d, can be given by the formula 



where xj/^ and F mn are the Fourier coefficients of \\r(x,y) and the (experimentally 
measurable) function F (x,y) = [M 2 I (Mx y My, R 2 ) - /J / (7 0 R'), respectively. 



In the general case of non-uniform intensity distribution on the object plane, the 
last two terms in equation (16) can be comparable in magnitude. However, in many real 
situations, where the illuminating wave has uniform intensity and absorption in the 
object is moderate, but not negligible, the last term in equation (16) may be smaller than 
25 the rest. In such cases equation (18) can be solved by perturbation methods. In 
particular, it is easy to show that an equation similar to (20), but with the constant 

intensity I 0 in the numerator of the function F (x,y) replaced by the measured variable 
intensity distribution I(x,y), can be considered as a first order approximation to the exact 
solution of equation (18). 



15 




(20) 



20 
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Finally, when the last term in equation (16) has a similar magnitude to the 
preceding one, system (18)-(19) can be solved numerically using one of the well known 
methods for the Laplace equation on the plane. 

5 

It is believed that the described derivation and radiographic method and apparatus 
provide a practical technique for deriving phase-contrast information contained in one or 
more detected intensity distributions in two dimensions. The configuration shown in 
Figure 1 is straight forward and akin to configurations traditionally used for absorption- 
10 contrast radiography. The disclosed method for deriving the phase-contrast image is 
particularly suited to planar detectors of a pixellated or grid construction or which are 
read in an incremental manner, and all the more suitable to the combination of such 
detectors and digital computers. 

15 Furthermore, by processing the detector values for absorption-contrast by 

conventional methods, one can simultaneously derive both phase-contrast and 
absorption-contrast images from, say, radiographic recordings in two or more two- 
dimensional planes D P D 2 . 

20 Recording of two image intensity distributions at different distances D r D 2 in 

order to carry out phase retrieval may not always be convenient or practical. For 
example, the measurement of intensity in two planes at different positions along the optic 
axis usually requires precision movements of the detector system D (indicated by arrow 
d in Figure 1) which may be undesirable in real-life applications. A second method of 

25 phase retrieval in the present approach may be achieved by measuring two or more 
intensity distributions in an image plane at fixed distance from the object but with two or 
more distinct distributions of incident energy, e.g. different wavelengths in a simple case 
or even more advantageously by using a 2-dimensional energy-sensitive detector such as 
those based on Cd-Mn-Te. With reference again to Figure 2, consider a paraxial 

30 wavefield in the free half-space z > 0 in an arbitrary state of temporal and spatial 
coherence with complex amplitude 
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t/(r,*)= Jf/(r,v)exp(/27rv0dv, (21) 

where v = ck/ {2k) is the frequency, k is the wavenumber and c is the speed of light in 
vacuum. Let S(r, v) be the (power) spectral density, then 

5 

5(r ± , /?, v) - S(r ± , v) = -/?V X • [5(r x , v)V 1 i^(r_ L , v)], (22) 

where /^(r^, v) = argf/(r ± , v) and we omit z=0 from the list of arguments of all 
functions. In the monochromatic case we have U(r, v) = U(r)5(v- v 0 ), 
10 S(r,v) = /(r)S(v- v 0 ) and the phase ky/(r,v 0 ) coincides with ky/(r ± ) = argC/(r ± ). 

The TIE can be obtained by integrating equation (22) over the frequency range: 
7(r ± , R) - /(r x ) — J?V A •[/ S(r ± , v) V x ^(r ± , v) d v] , (23) 

15 

where 7(r) = j S(r. v)dv is the (time-averaged) intensity [14]. As already discussed, this 

form of the TIE relates the difference between the intensity distributions in the object 
and the image planes to the phase derivatives in the object plane. 

20 If spatial variations of the spectral density are much weaker than the 

corresponding phase variations, then equation (22) can be simplified: 

5(r 1 ,/?,v) = 5(r J _,v)[l-/J V>(r x ,v)]. (24) 

25 This situation is typical for phase objects, i.e. objects with the real part of the 

refractive index varying much stronger than its imaginary part at the given wavelength. 
This is the case e.g. for thin biological samples imaged with high-energy x-rays. For 
such objects conventional images based on absorption often have poor contrast, while 
phase-contrast images are much more informative [5, 8, 9]. 
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Let us assume that the propagation of x-rays through the sample is paraxial and 
rectilinear. Then the spectral density and the phase in the object plane can be expressed 
as 

5 

S(r ± , v) = S in (r ± , v)exp {-[ ju(r ± ,z\ v)dz'} , (25) 

y/(r ± , v) = y/ in (r ± , v) - J 5(r ± ,z v)dz ' , (26) 

10 where the values with superscript "in" correspond to the wave incident on the sample, 
integration is through the sample along the lines parallel to the optic axis, 1 - 8 is the real 
part of the refractive index and fi is the linear absorption coefficient. If all wavelengths 
present in the spectrum are sufficiently far away from the absorption edges of the sample 
material, then 

15 

)U(r,v) = <T 3 //(r,v 0 ), (27) 

<5(r,v)-cr 2 (5(r,v 0 ). (28) 

20 where a = v 0 / v and v 0 is some fixed frequency from the spectrum. Substituting 
equations (25)-(28) into equation (24), we obtain 

S(R,v) = S in (v)exp(-<J 3 M 0 )[l- R V^ /w (v) + /?<r 2 V^y/J, (29) 

25 where we omitted r ± from the list of arguments of all functions and denoted 

M 0 = j fJ.(r ± ,z',v 0 )dz' and y/ 0 = jS(r ± ,z\ v 0 )dz' . The last two terms in the square 

brackets in equation (29) must be much less than one if for example the paraxial 
approximation is valid. 
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If the characteristics of the radiation incident on the sample, i.e. S m (v) and 
y/ in (v), are known, and the spectral density at a given frequency at the image plane, 
S(R,v), can be measured (e.g. with the help of an energy-sensitive detector), then 
equation (29) can be used for phase retrieval from two such measurements at different 
: wavelengths. Indeed, if, for example, S(R,v) and S(R,v 0 ) have been measured, then 
writing e.g. equation (29) for v and v 0 and excluding exp{-cr 3 M 0 } from these 
equations, it is possible to derive 



v'u, i-y-/?[Vlv/ f "(v)-yc7 3 vy(v 0 )] 



10 



where 



y= S(Jtv) 



S in (v 0 ) 
S in (v) |_S(/?,v 0 ) 



(31) 



15 The phase y/ 0 (r L ) can be recovered by solving the Poisson equation (30). Note 

that any uniform boundary conditions can be assigned to y/ 0 , if the image is surrounded 
by the unperturbed primary beam. The solution to equation (30) with typical uniform 
boundary conditions is always unique at least up to an arbitrary additive constant. It is 
obvious from equations (30) and (31) that the stability of phase recovery is determined 

20 by the ratio a = v 0 / v of the two frequencies and the ratio of the spectral densities in 
equation (31). In an experiment these values can be chosen to be sufficiently different 
from unity, in which case the recovery of the phase using equation (30) will be stable 
with respect to errors in the experimental data. Therefore, this technique may have an 
advantage over the earlier-described method of phase recovery from intensity 

25 distribution in the beam cross-sections at two different positions along the optic axis. 
Indeed, in this latter case the z-distance between the cross-sections must be small and the 
difference between the two intensity distributions has to be small too, so the calculation 
of the phase Laplacian is a numerically unstable procedure of the "division of zero by 
zero" type. 
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As with the earlier embodiments of the invention, an analog of equation (30) 
suitable for in-line phase imaging with a small polychromatic source can be derived. The 
additional curvature appearing in the wavefront due to the quasi-spherical nature of the 
5 incident beam can be explicitly taken into account. Consider the following "quasi- 
spherical" analog of equation (22): 

M 2 S{Mv L ,R, v) - S(r ± ,v) = -*'V ± • [5(r ± , v)V x ^(r ± , v)] }, (32) 

10 where M = (/?, + R) I R { is the magnification coefficient, R' = (R~ l 4- R~ l ) _1 and R { is the 

source-to-object distance. Starting from equation (32) and proceeding exactly as above, 
it is straightforward to derive the Poisson equation for the phase similar to equation (30), 
but with y/ in (r x ,v) replaced by i// , "(r ± ,v) = i/A /w (r x ,v)-r^/(2/? 1 ), R replaced by R' 
and 7 replaced by 



15 



Y > = M w-*>) S(Mr ±J R,v) 



S in (r ± ,v 0 ) 



(33) 



In a slightly different experimental arrangement without an energy-sensitive 
detector it may be possible to perform intensity measurements in the image plane with 
20 two different monochromatic incident beams with frequencies v and v 0 . Then the above 
equations (30) and (32) still hold with the replacement of all spectral densities by the 
corresponding intensities. 

It is now proposed to outline a method for phase retrieval from polychromatic 
25 images obtained using only a conventional (not energy-sensitive) detector. Let us assume 
that the transversal variations of absorption are so weak that 

exp { -cr 3 M 0 (r ± ) } s exp { -a 3 M 0 } [1 - a" M' 0 (r L )] , (34) 
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where M 0 is the "average" absorption in the sample at v = v 0 , while \cr 3 M' 0 (r L )\« 1 
everywhere. The constant M Q can be determined from the knowledge of total intensity 
incident on the sample and the total intensity of the image. The weak local variations 
M 0 '(r x ) of absorption are assumed unknown. 

In the first crude approximation one may choose to neglect the term A/ 0 '(r ± ). 
Then the following Poisson equation for the phase can be easily derived from equation 
(29): 



10 -V> 0 = J - ? , (35) 

1 R\s out {v)o 2 <\v 



using the notation S out (v) = S in (v)zxp{-<j 3 M 0 ) . Therefore, in this case the phase can be 
found from only one polychromatic image I(R) = I(r ± ,R) = J 5(r x ,-/?, v)dv by solving 

equation (35), providing the properties of the incident radiation, i.e. S in (v) and y/ in (v), 
15 are known a priori. 



It is now proposed to demonstrate that, if the weak local variations M 0 '(r ± ) of 
absorption have to be taken into account, the phase can be found from two 
polychromatic images obtained with the incident beams having different spectral 
20 composition. Substituting equation (34) into equation (29) and integrating over the 
frequency range, we obtain 

a - bM; - cV^Vo + dM 'oV 2 ±Vo = 0, (36) 

25 where a = I(R)~ J>"'(v)[l - RV 2 ± y/ in (v)]dv 9 b = -J> w '(v)[l - /?V>'"(v)]<J 3 dv , 

c = Rjs OUI (v)cr 2 dv and d = -j S out (v)c 5 dv . Given two images, /,(/?), 7=1,2, measured 
with incident beams having different spectral composition (described by the a priori 
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known distributions S'"(v) and y^"(v), 7=1,2), M 0 ' can be eliminated from the 
corresponding pair of equation (36) and obtain 



_v> 0 = eA^A , (37) 

b x c 2 — b 2 c x 4- a 2 (d y — d 2 b x I b 2 ) 



with all the quantities defined for each of the incident beams exactly as above, e.g. 
cij =Ij(R)- j S° ut (v)[l - RVly/t" (v)]dv , j=l,2, etc, similarly, it is possible to solve for 
M o to obtain an effective pure absorption-contrast image. 

10 These expressions have demonstrated that the x-ray phase can be recovered using 

the Transport of Intensity equation from polychromatic images obtained at a fixed 
position along the optic axis. It can be done with a quasi-plane or a quasi-spherical 
paraxial incident wave, whose spectral density and the phases of the monochromatic 
components are known a priori. The approach would preferably involve a substantially 

15 complete characterisation of the source prior to phase retrieval experiments. The above- 
described embodiments of the invention demonstrate that effective phase- and 
absorption-contrast images can be obtained, in a variety of ways including any one of the 
following: 

(1) From two images obtained with monochromatic incident beams with two 
20 different wavelengths (photon energies). 

(2) With a polychromatic incident beam and measurement of the image 
intensity as a function of x-ray energy using an energy-sensitive detector (e.g. 
based on cadmium manganese telluride), with some binning of the energy values 
to give sufficient image intensity in at least two energy ranges (bins). 

25 (3) With a polychromatic incident beam and measurement of the image 

intensity using some energy selectivity in the detector, for example by the use of 
energy filters based on the use of absorber foils. In practice this might involve 
using a pair of x-ray films with an appropriate filter foil between them to 
simultaneously record suitable image intensity data. 
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(4) From two polychromatic images obtained with incident beams having 
different spectral composition and a conventional, rather than energy-sensitive, 
detector. 

5 In each of these above situations, just one measurement of an image is sufficient 

for phase recovery if the transverse spatial variations of absorption are negligible. 
Otherwise, it is possible to obtain separate effective phase- and absorption-contrast 
images from the image intensity data. 

10 Figure 3 diagrammatically illustrates a variety of arrangements for deriving plural 

intensity records for respective different energy distributions. Such measurements can 
be performed, for example, with an energy-sensitive detector D E , or with a source 
tunable by adjustable control C to vary its wavelength profile, or with a polychromatic 
source S and interchangeable monochromators or filters either ahead of the object (F,) or 

15 after the object, e.g. as a shutter F 2 front of the non-energy sensitive CCD detector D. If 
absorption in the sample is negligible, the phase can be found from intensity of only one 
polychromatic image, provided the properties of the incident radiation are known a 
priori. 

20 Other approaches to using two or more different spectral distributions (i.e. 

alternatives to approaches based on transport-of-intensity equations) may involve 
deconvolution of the spectral distribution from the data or, alternatively, a least squares- 
type fitting procedure to extract Vfy. Even more sophisticated methods based on 
Bayesian analysis might also be used. It is proposed to now briefly outline phase and 

25 absorption retrieval by a technique based on a Fourier optics approach. 

The treatment of phase retrieval in the spherical-wave case is analogous to the 
treatment for the plane-wave case and quite straightforward. This section is thus 
confined to the latter. 

30 
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It may be shown that in the plane-wave case, Cowley's form of the Kirchhoff's 
formula for the x-ray wave function can, with certain small-angle approximations, be 
reduced to: 

^(x,y) - * 0 < x ,y) exp(-iki^) * q(x,y) (38) 

5 

and therefore 

W(x,y)] -*(f..f y ) = * 0 (*.y) expti^rX^f 2 ^f y 2 )] Q(f x ,f y ), (39) 

10 where we adopt the convention that real-space functions are in lower-case and the 
corresponding frequency-space functions (after taking Fourier transforms) are in upper- 
case, e.g. Q(f f ) = ^[q(x,y)]. For a further discussion of the steps involved in this 
derivation, see J.M. Cowley (1975) Diffraction Physics (North-Holland: Amsterdam.). 
Relevant theoretical discussion is also to be found at Pogany et al (1997) Rev. Sci. Inst. 

15 68,2774-2782. 

The transmission function can be written, assuming <|>(x,y)t(x,y) and |x(x,y)t(x,y) 
are sufficiently small, as q(x,y) ~ 1 + i<|) , (x,y)-fi , (x,y), where <|>Xx,y)=-<t)(x,y)t(x,y) and 
|j/(x,y) = (i(x,y)t(x,y)/2. The Fourier transform of the transmission function is then given 
20 by Q(f x ,f ) ~ 8(f ,f) + iOXf ,f y )-M'(f ,f). Equation (39) can then be written as 



*(f x ,f y ) - ^ exp[ix(f x ,f y )] [6(f x ,f y )-i*'(f x ,f y )-M'(f x ,f y )l, (40) 

25 where %(f ,f ) = nXR 2 (f x 2 +f 2 ). By expanding (40) and taking inverse Fourier transforms 
of both sides, one can obtain, to first order in $ and |a\ and omitting the functional 
dependence: 

I — 1 - 2<j>' * ^- l [sin( x )] - 2p' * ^"'[cosCx)]. (41 ) 
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Consequently, in the case of a pure-phase object, one can retrieve the <|>t- 
distribution from a single image: 



and, in the case of a pure-absorption object: 



10 



In the more usual case of an object for which both phase and absorption effects 
15 are significant, two images, I, and I 2 , which have been collected with different values of 
R 2 and/or x-ray energy (we are discussing the monochromatic case here, but in the 
polychromatic case two different voltage settings of the x-ray source, giving two 
different spectral distributions, would be appropriate) could be used. Assuming, in the 
absence of any absorption edges for constituent elements, that \i 2 = (kJX^^ and that ty 2 
20 = (KJX^ty^ the resulting simultaneous equations can be solved to yield: 





5T(I 2 


-l)cos(x,) - 


V 
\ 


3 

^-(I,-l)cos(x 2 ) 


2 


K 
\ 

— V J 


3 

sin(xj) cos(x 2 ) - 


K 

-2 sin(x 2 ) cos(xi) 



(44) 



and 
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fit = & A 





V 
\ 


^-(1,-1) sin(x 2 ) - ^-(I 2 -D sin( Xl ) 


V 


3 

sin(x,) cos(x 2 ) - 


V 
\ 


sin(x 2 ) cos(x,) 



5 The advantages of using the multiple energy method of retrieving phase rather 

than the multiple distance method include the following: 

• ability to record data by electronic rather than mechanical change of 
imaging conditions (for example by rapid switching of tube voltage or 
energy-sensitive detectors) leading to the possibility of very rapid phase 

10 retrieval; and 

• ability to use 2-dimensional energy dispersive detectors (e.g. those based on 
CdMnTe) as these become available, in order to record a large number of 2- 
dimensional images for different x-ray energy ranges. 

15 Example 1 

To test the validity of the method of the invention utilising transport-of-intensity 
equations, a numerical example was prepared with simulated data, and the successive 
steps are illustrated in Figures 4 to 9. 

20 

A perfectly monochromatic x-ray point source was assumed with wavelength X = 
0.154 nm (CuK a radiation), and source-to-object and object-to-image distances R { = R 2 = 
0.2m, so that R' = 0.1m and M = 2 (Figure 2). A pure phase object was also assumed, of 
size 640 x 640 \im 2 , producing the phase distribution on the object plane (x,y f 0) shown in 
25 Figure 4. The range of the phase values was [-0.8667, 0.6706] radians. 

The intensity distribution I(x,y,R 2 ) on the image plane was then calculated by 
computing the Fresnel integral (7 above) on a grid with 128 x 128 pixels. The scaled 
image M 2 I(x / M,y / M, R 2 ) of the obtained intensity distribution 7(x,^,R 2 ) is presented in 
30 Figure 5. For comparison the calculated distribution of the phase Laplacian, -V 2 (p(x,y,0) 
was plotted and is shown in Figure 6. Obviously, the scaled distribution of intensity on 
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the image plane and the distribution of the Laplacian of the phase on the object plane 
appear almost identical, which supports formula (16) in the case of negligible absorption. 

At the next stage the computed intensity distribution on the image plane and 
5 formula (20) were used to reconstruct the phase. The result of this reconstruction is 
shown in Figure 7. Again, it is in a very good agreement with the original phase from 
Figure 4, with the relative mean-square error of this reconstruction equal to 4.65%. 

The performance of the method was also tested in the case of non-uniform 
10 intensity, i.e. when the last term in formula (16) was significant. An intensity 
distribution on the object plane was simulated according to the formula: 
I( Xf y) = -1.647 + 3 x exp {-(x 2 +y ) / (2*1280 2 )}. 

The values of the intensity had the range of [1, 1.1243]. The intensity 
15 distribution /(Ot,y,R 2 ) on the image plane was computed by evaluating the Fresnel 
integrals. The non-uniformity of intensity distribution resulted in the "vignetting" of the 
corresponding image (Figure 8), which clearly reflects the non-uniformity of intensity on 
the object plane. However, by subtracting the intensity distribution on the object plane 
from that on the image plane, dividing the difference by the scaled interplanar distance 
20 R' = 0.1m and applying formula (2) we were able to reconstruct the original phase with a 
root mean-square error of 6.67% (Figure 9). This result confirms that in the case of 
moderate non-uniformity of intensity distribution on the object plane (e.g. weak 
absorption effects and/or weak non-uniformity of the incident beam), the described 
simple variation of formula (20) can be effectively used for quantitative retrieval of the 
25 phase shifts produced by the object. 

Example 2 

Figures 10 to 15 show an example of phase and absorption retrieval from two 
30 images at different values of R 2 . The CSIRO logo has been used as the absorption (jit) 
and phase (c(>t) distributions (Figures 10, 11 respectively). For the former the values are 
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either 0.0 (black) or 0.1 (white), and for the latter (inverted logo) -0.1 (black) or 0.0 
(white). The images (intensity distributions) were calculated for the plane-wave case, 
lA radiation, 0.2x0.2 mm 2 (512x512 pixels) area, and R 2 = 80cm (Figure 12) and R 2 = 
160cm (Figure 13). These conditions are such that the images are not immediately 
5 interpretable. Applying the retrieval algorithm (44,45) discussed, we obtain the results 
shown (Figures 14, 15) which while not perfect, are clearly recognisable. 
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CLAIMS: 

1. A method of obtaining an image of the phase change introduced by an object in 
penetrating radiation incident on the object, including: 

5 irradiating the object with penetrating radiation having high lateral spatial 

coherence; 

receiving at least a portion of said radiation at detector means after the 
radiation has emerged from the object and thereby obtaining and storing at least 
two intensity records for the received radiation each including intensity values at 
10 predetermined intervals; and 

utilising these values to derive a grid of values defining an image of the 
phase change introduced by the object in the penetrating radiation; 

wherein said intensity records are obtained at a uniform finite distance 
after the radiation has emerged from the object, and are for respective different 
15 energy distributions of the detected radiation. 

2. A method according to claim 1 wherein the respective different energy 
distributions are obtained by altering the energy spectrum of the radiation irradiating the 
object. 

20 

3. A method according to claim 1 wherein the respective different energy 
distributions are obtained by providing for said detector means to provide intensity as a 
function of energy in a certain energy band or band(s). 

25 4. A method according to claim 1, 2 or 3 wherein said derivation includes solving 
one or more differential transport-of-intensity equations relating the phase at a plane of 
the object to the evolution of the intensity distribution along the direction of propagation, 
utilising predetermined uniform boundary conditions. 

30 5. A method according to claim 1, 2 or 3 wherein said derivation includes solving 
Fourier-optics equations. 
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6. A method according to any preceding claim wherein said intensity values also 
reflect absorption contrast in the object, and the method further includes utilising said 
values to derive a grid of values defining an effective pure absorption-contrast image of 
the object. 

5 

7. A method according to any preceding claim wherein said penetrating radiation 
comprises x-ray radiation. 

8. A method according to claim 7 wherein said x-ray radiation is in the range 0.5 
10 KeV to 1 MeV. 

9. A method according to claim 7 or 8 wherein said irradiating radiation is 
substantially monochromatic. 

15 10. A method according to claim 7 or 8 wherein said irradiating radiation is 
polychromatic. 

11. A method according to any one of claims 7 to 10 wherein said irradiating 
radiation is from a substantially point source of full width half-maximum 40 |iim or less. 

20 

12. Apparatus for obtaining an image of the phase change introduced by an object in 
penetrating radiation incident on the object, including: 

means to provide a source for irradiating an object with penetrating 
radiation having high lateral spatial coherence; and 
25 detector means for receiving at least a portion of said radiation after the 

radiation has emerged from the object whereby to generate at least two intensity 
records for the received radiation each including intensity values at 
predetermined intervals; 

wherein said detector means is arranged for obtaining said intensity 
30 records at a uniform finite distance after the radiation has emerged from the 
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object, and energy characterising means is provided whereby said intensity 
records are for respective different energy distributions of the detected radiation. 

13. Apparatus according to claim 12 wherein said energy characterising means is 
5 arranged to alter the energy spectrum of the radiation irradiating the object 

14. Apparatus according to claim 13 wherein said energy characterising means 
includes means rendering the detector means able to provide intensity as a function of 
energy in a certain energy band or band(s). 

10 

15. Apparatus according to claim 12, 13 or 14 further including a computer program 
product having a set of machine readable instructions which, when installed in a 
computer having a suitable operating system and memory means, configures the 
computer to be operable to utilise said values to derive a grid of values defining an image 

15 of the phase change introduced by the object in the penetrating radiation. 

16. Apparatus according to claim 15 wherein said derivation includes solving one or 
more differential transport-of-intensity equations relating the phase at a plane of the 
object to the evolution of the intensity distribution along the direction of propagation, 

20 utilising predetermined uniform boundary conditions. 

17. Apparatus according to claim 15 wherein said derivation includes solving 
Fourier- optics equations. 

25 18. Apparatus according to any one of claims 12 to 17 further including a source of 
x-ray radiation as said source for irradiating the object. 

19. Apparatus according to claim 18 wherein said x-ray radiation is in the range 0.5 
KeV to 1 MeV. 

30 

20. Apparatus according to claim 18 or 19 wherein said irradiating radiation is 
substantially monochromatic. 
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21. Apparatus according to claim 18 or 19 wherein said irradiating radiation is 
polychromatic. 

5 22. Apparatus according to any one of claims 18 to 21 wherein said source is a 
substantially point source of full width half-maximum 40 |Lim or less. 

23. A method of obtaining an image of the phase change introduced by an object in 
penetrating radiation incident on the object, from one or more two-dimensional intensity 

10 records of penetrating radiation after it has traversed the object, the radiation being of 
high lateral spatial coherence when incident on the object and the or each record being 
obtained at a finite distance after the radiation has emerged from the object incorporating 
phase-perturbed components within a surrounding field of radiation either uniformly 
phase-perturbed or not phase-perturbed, the method including: 

15 storing intensity values from the or each record, at predetermined 

intervals; 

utilising these values and any predetermined uniform boundary conditions 
to derive a grid of values defining an image of the phase change introduced by 
the object in the penetrating radiation, by solving a differential transport-of- 
20 intensity equation relating the phase at a plane at the object to the evolution of the 

intensity distribution along the direction of propagation. 

24. A method of obtaining an image of the phase change introduced by an object in 
penetrating radiation incident on the object, including: 

25 irradiating the object with penetrating radiation having high lateral spatial 

coherence; 

receiving at least a portion of the said radiation at a detector at one or 
more finite distances after the radiation has emerged from the object 
incorporating phase-perturbed components within a surrounding field of radiation 
30 either uniformly phase-perturbed or not phase-perturbed, and thereby obtaining 
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and storing intensity values for the received radiation at predetermined intervals; 
and 

utilising these values and any predetermined uniform boundary conditions 
to derive a grid of values defining an image of the phase change introduced by 
5 the object in the penetrating radiation, by solving a differential transport-of- 

intensity equation relating the phase at a plane at the object to the evolution of the 
intensity distribution along the direction of propagation. 

25. A method according claim 23 or 24 wherein said intensity values also reflect 
10 absorption contrast in the object, and the method further includes utilising said values to 

derive a grid of values defining an effective pure absorption-contrast image of the object. 

26. A method according to claim 23, 24 or 25 wherein said penetrating radiation 
comprises x-ray radiation. 

15 

27. A method according to claim 26 wherein said x-ray radiation is in the range 0.5 
keV to 1 MeV. 

28. A method according to claim 26 or 27 wherein said irradiating radiation is 
20 substantially monochromatic. 

29. A method according to claim 26 or 27 wherein said irradiating radiation is 
polychromatic. 

25 30. A method according to claim 29 wherein said equation includes a spectrally 
weighted term or factor dependent on the square of the respective wavelength 
components. 

31. A method according to any one of claims 23 to 30 wherein said boundary 
30 conditions include uniform Dirichlet, Neumann or periodic boundary conditions, and are 
selected to achieve a unique solution of the equation for phase, at least up to an arbitrary 
constant component. 
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32. A method according to claim 31 wherein the solution further utilises one or more 
optical conditions selected from the group consisting of a small wavefront curvature for 
the incident radiation, absence of focal points between the object and image, and uniform 

5 illumination of the object. 

33. Apparatus for obtaining an image of the phase change introduced by an object in 
penetrating radiation incident on the object, including: 

a source for irradiating an object with penetrating radiation having high 
10 lateral spatial coherence; 

a detector for receiving at least a portion of said radiation a finite distance 
after the radiation has emerged from the object incorporating phase-perturbed 
components within a surrounding field of radiation not phase-perturbed or 
uniformly phase-perturbed , whereby to generate intensity values for the received 
15 radiation at pre-determined intervals; and 

computer means, including a stored program of machine readable 
instructions, operable to utilise said values and any predetermined uniform 
boundary conditions to derive a grid of values defining an image of the phase 
change introduced by the object in the penetrating radiation, by solving a 
20 differential transport-of-intensity equation relating the phase at a plane at the 

object to the evolution of the intensity distribution along the direction of 
propagation. 
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